Centralize the equation-of-state expressions in the Riemann solvers - #1762

Gives each equation-of-state expression one definition, so that adding a second EOS (#1638) is a
change inside one module rather than an edit to every solver.

The operators

s_compute_mixture_coefficients(alpha_rho, alpha, rho, gamma, pi_inf, qv)
s_compute_energy(pres, alpha_rho, alpha, vel_sum, E)
s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)   ! _avg for interface averages
f_bulk_modulus(pres, gamma, pi_inf)
f_pressure(e_int, gamma, pi_inf, qv)
f_phase_internal_energy(pres, alpha, alpha_rho, gamma, pi_inf, qv)
f_isentrope_exponent(gamma) / f_isentrope_pressure(pi_inf, gamma)

s_compute_speed_of_sound used to take H, |u|^2 and qv and then undo them; for a real state all
three cancel, so it takes none of them and #1707 becomes unrepresentable.

The f_ operators take coefficients rather than a fluid index, so a mixture - whose effective gamma
and pi_inf already come from s_compute_mixture_coefficients - is the same call as a single fluid.
That is what lets the sound speed use the bulk modulus: every one of its cases is a bulk modulus over a
density, differing only in how the phases mix (Wood's law is their harmonic mean, the six-equation model
their arithmetic mean, the rest the mixture value). Written that way the mpp_lim sub-branch of the
bubbles case is visibly the same expression as the plain case and collapses. The four rewrites were
checked against the originals over 200k random states, agreeing to three ulp.

gs_min and ps_inf already are the Tait exponent and reference pressure, 1/gamma + 1 and
pi_inf/(1 + gamma). Three places rebuilt them by hand - m_qbmm's pi_infs(1)*(n - 1)/n reduces
exactly to ps_inf(1) - and HLLC's star state now reads as the relation it is, p* = (p + B) xi**n - B.

Kinetic, magnetic and elastic energy stay at the call sites; they are not EOS terms.
s_compute_pressure decides what the internal energy is, then inverts once.

Copies collapsed: mixture rule 9, bulk modulus 7, pressure inversion 5, elastic strain energy 4,
per-phase internal energy 3, isentrope parameters 3.

A second sweep for inline stiffened-gas arithmetic found four more, all now routed through the
operators: post_process's three derived quantities (specific heat ratio, liquid stiffness, sound
speed - the last sat beside a call already using f_bulk_modulus); s_compute_cson_from_pinf, a
fifth sound-speed copy in the E/H round-trip form #1707 was about, which built E and H only
to cancel |u|^2 - a real cancellation loss, not just verbosity, and its deletion also removes a
USING_AMD workaround and two arguments from a [seq] device routine; and the pressure-based wave
speed in HLL and HLLC, whose (0.5 + Gamma)/(1 + Gamma) is exactly Toro's two-shock coefficient
(n + 1)/2n and whose (p*/p - 1)*p divided by p only to multiply it back.

Behaviour changes

Each is a forward/inverse disagreement or a wrong mixture, none exercised by any test:

• hll and lf lacked the bubbles_euler special case, computing gamma = alf*gammas(1) with
  alf a void fraction. They now take HLLC's rule, which is the published one: under bubbles_euler
  the last advection slot is the void, not a material (eqn_idx%alf == eqn_idx%adv%end), so the
  mixture sums the material slots only and a single carrier liquid keeps its own coefficients
  undiluted - Gamma_l p_l = (E - rho|u|^2/2)/(1 - alpha) - Pi_inf_l, the void entering only through
  the 1/(1 - alpha) that s_compute_pressure already applies. HLLC and m_acoustic_src are
  unchanged by this; a two-fluid bubbly case is added, since nothing ran bubbles_euler above one
  fluid. num_fluids > 2 is now refused at validation: the closure is written for one carrier liquid, the
  two arms above it disagreed, and nothing exercised either (bubbles_euler with num_fluids > 2 runs an unvalidated mixture closure #1786).
• m_cbc had the same gap in derivative form. It differentiates the mixture rule with no
  bubbles_euler case, so three of its four coefficient derivatives were wrong there - they should
  vanish, because the rule makes gamma, pi_inf and qv constants.
• HLLC's bubbles energy dropped qv while s_compute_pressure subtracts it.
• m_ibm dropped qv twice (ghost-point energy and per-phase internal energy), put the dynamic
  pressure inside the (1 - alpha) scaling instead of outside, and used alpha_IP(1) where
  eqn_idx%alf is alpha_IP(num_fluids).
• post_process's energy budget formed sum(alpha*alpha*rho) rather than sum(alpha*rho); that
  density feeds the sound speed, so its Mach number was wrong wherever a volume fraction was below one.
  Its gas internal energy Egint was alpha_g*gamma_g*p, missing the alpha_g*pi_inf_g + alpha_rho_g*qv_g
  that make up the rest of a stiffened-gas internal energy; it now calls f_phase_internal_energy, so it
  agrees with the pressure inversion that produced p. Identical for a gas phase with pi_inf = qv = 0.
• The hypoelastic probe pressure rebuilt the kinetic energy from one momentum component while the
  correct dyn_p was passed in and ignored - wrong in 2D and 3D. mom is gone from
  s_compute_pressure, so no caller can reintroduce it.
• m_acoustic_src never applied the subgrid void correction. It was the last site hand-rolling a
  sound speed. The Tait form it used is the same stiffened-gas expression as f_bulk_modulus - with
  n = 1/Gamma + 1, (n-1)/n = 1/(1+Gamma) - so that half is pure de-duplication, and the downstream
  c/(small_gamma - 1) was c*Gamma written the long way round. The behaviour change is the dilution:
  a wave injected into a bubbly medium travels at the mixture speed c_l/(1 - alf), and the hand-rolled
  code used the undiluted carrier-liquid speed. 39 goldens move, all bubbles_euler, all by ~1.7% at
  the source cell - that is 1/(1 - alf) at alf = 0.04, and it is visible only because the test grids
  are deliberately coarse (49x39 in 2D against 299 in 1D), so a 2% shift at the source sits well outside a
  1e-10 tolerance. Lagrange and non-bubbly acoustic cases do not move: the guard excludes them.
  The correction is inlined rather than calling s_compute_speed_of_sound, which faults CCE OpenACC from
  that loop for reasons not yet isolated (CCE OpenACC faults when the acoustic loop calls s_compute_speed_of_sound #1794); the comment there says to fold it back once they are.
• The subgrid void fraction is now named where it is used (Bubble sound speed: alpha -> 1 should assert, not silently produce Inf #1793). The bubble branch of the mixture
  sound speed is c = c_l/(1 - alf) - the carrier-liquid speed with an O(alf) correction, the same c_l
  that m_qbmm forms for Keller-Miksis. alf is dilute by construction, so a value near one is a wrong
  index or a case outside the expansion's regime, not a number to clamp; the toolchain warns at case load
  instead. No arithmetic changes.
• relativity silently ignored pi_inf. The relativistic enthalpy h = 1 + (Gamma + 1)p/rho
  carries no stiffness, and it was written out twice - once in m_riemann_solver_hll under a
  ! Hard-coded EOS comment, once inside s_compute_speed_of_sound's own relativity branch. It is
  now f_relativistic_enthalpy, and a stiffened fluid with relativity is refused at validation
  rather than quietly given an ideal-gas answer.

Test comparison

A golden that records an analytically zero field stores only roundoff, and comparing it pointwise tests
the compiler's association order: 421A6AD9 failed on a sign flip of a 5.8e-07 crumb. A field whose
every golden value already lies within the absolute tolerance of zero is now checked for staying in a
band rather than for reproducing its residue. This also closes a hole - compute_error returns NaN
when the golden value is exactly 0 and is_close passed NaN unconditionally, so ~2800 fields were
not checked at all.

The band is the larger of the tolerance and ten times the golden's own magnitude, not the tolerance
alone. Held to the bare tolerance the rule misfires on a field that merely sits near it: single
precision scales QBMM's 1e-10 by 1e8, so its tolerance is 1e-2 and its bubble fields are
constants at 0.9999997 of that - read as zero, then failed for drifting a few ulps past it. Scaling
by the golden judges such a field against itself. A field storing genuine roundoff is orders below the
tolerance, so its band is unchanged: the 5.8e-07 sign flip still passes and a candidate growing from
roundoff to 1e-3 still fails.

Verification

MI210, AFAR amdflang, OpenMP offload. 5eq_rk3_weno3_hllc, three runs each, one node per comparison:

0 regression(s) across 471 kernels in static GPU resources, so the [seq] operators inline even in
the per-phase, per-face HLLC flux loop.

Full suite 684 passed. Targeted set of 196 covering every riemann_solver, avg_state,
wave_speeds, low_Mach and mixture_err combination plus chemistry, reactive burn, phase change
and hypoelasticity - 196 passed, 0 failed. The only golden regenerated is the new 2D probe case,
for #1773 below.

The nvfortran fix was verified as a before/after pair on one machine, nvfortran 24.1 in a container:
the pre-fix tree reproduces fort2 TERMINATED by signal 11 on the HLL solver and the post-fix tree
builds clean under both --gpu mp and --gpu acc. Reverting only the molecular_weights change
brings the crash back, which is what attributes it. The two CCE fixes were confirmed on Frontier -
the device-global rule by a six-build bisection, and the runtime fault by all four CCE lanes passing
once the thermochem calls moved to the loop body.

Scope

Deliberately left. Each says what it would take, so none of it needs rediscovering:

• m_phase_change's caloric relations. Entropy, enthalpy and e(p,T) appear once each inside
  s_infinite_ptg_relaxation_k. Extracting them one at a time would be net-negative today, and a
  mechanical-only EOS supplies none of them, so the first question is a validator prohibition rather
  than code. Filed as Phase change assumes a caloric equation of state that a mechanical EOS does not supply #1784.
• m_pressure_relaxation's isentrope inversion. Needs rho(p) and drho/dp in closed form;
  JWL's isentrope has neither, so it needs a numerical inverse inside a per-cell Newton loop. The one
  genuinely new piece of solver code multi-EOS requires. Filed as Pressure relaxation inverts the stiffened-gas isentrope in closed form; a generic EOS needs a numerical inverse #1785.

Supersedes #1714 (same goal; its type(eos_state) cost 20% because a derived-type dummy on a [seq]
device routine forces a memory ABI). Closes #1707, #1708, #1715, #1769, #1773, #1774, #1777, #1778, #1781, #1786. Unblocks #1638.

Also filed from this work, not fixed here: #1784, #1785, #1793, #1794, #1779 (fixed upstream by #1780),
and the audit comments on #1682, #1687 and #1693 correcting their stated impact.

m_acoustic_src was the last site hand-rolling a sound speed; it now uses f_bulk_modulus and applies
the subgrid void correction. See the behaviour-change entry above. The remaining open item is the CCE
mechanism (#1794), not the physics.

Bisect caveat: a94f914a is not independently runnable. Its case_validator.py calls
check_eos_selector, whose CONSTRAINTS["fluid_pp(1)%eos"] entry does not arrive until f009397c,
so ./mfc.sh at exactly that commit raises KeyError. The commits were split by file rather than by
hunk. The branch tip is correct; only a bisect landing on that one commit sees it.

Folded in: the Riemann macro cleanup (#1769)

The inline_riemann.fpp / inline_capillary.fpp macros inline computational statements into the
solver loops and bind to the caller's locals by name. A missing local is reported at the call site
rather than in the macro, the bodies do not appear when reading the solver, and every variable they
touch has to be carried by hand in the GPU private() lists.

All nine are now procedures, and both include files are deleted:

inline_riemann.fpp and inline_capillary.fpp are both gone. The Lax-Friedrichs solver expanded
compute_average_state and used none of it, so that call site is deleted rather than converted -
which is what #1715 was reporting.

The low-Mach one is the interesting case. It branched on riemann_solver and on low_Mach at
runtime, but every one of its eight call sites already knew both answers: low_Mach = 2 requires
riemann_solver = 2, HLLD forbids low_Mach entirely, and each site sits inside its own low_Mach
guard in a file that is one solver. Every branch in the macro was dead at every site, so it was split
by leaf. The two pcorr functions absorb the low_Mach == 1 test and return zero otherwise, which
drops the repeated else pcorr = 0 from five sites; the low_Mach == 2 guard stays visible at its
three sites because that leaf mutates the wave-normal velocities before the wave speeds are
computed. zcoef, vel_L_tmp and vel_R_tmp become local and leave five GPU private lists.

Static GPU resources: 0 regression(s) across 471 kernels for the low-Mach change, measured as a
before/after pair on one node.

Compiler portability

Extracting solver code into $:GPU_ROUTINE device routines is not compiler-neutral. Two constructs
built cleanly under amdflang OpenMP offload and CCE OpenMP, and broke elsewhere:

CCE OpenACC rejected the pressure-relaxation solver:

ftn-7066 ftn: ERROR S_CORRECT_INTERNAL_ENERGIES, m_pressure_relaxation.fpp, Line = 211
  Unsupported accelerator code error: Global in accelerator routine without declare -- num_fluids

Six builds on Frontier pinned the rule down. An array whose bound is a device global may be passed
to a device routine from a parallel-loop body, but not from inside another acc routine seq.
Three plausible readings are wrong, each measured rather than argued: it is not call depth (moving
the construct up one level only moved the error), not the optional dummies (an operator with none
still failed), and not the local's bound (declaring it dimension(num_fluids_max), a parameter,
still failed). A variant keeping the num_fluids locals, the gather loop and every in-module call
built once the one cross-routine array pass was removed.

The reported line is where the compiler gave up, not the cause: remove one trigger and the message
walks forward to the next call. And only the plain lanes fail - under --case-optimization those
bounds are parameters, so the green Case Opt lane beside a failing plain one is the signature.

Every accepted call site in the tree already obeys the rule - m_cbc, m_ibm, m_bubbles_EL and
s_compute_cell_state all call this kernel from a loop body - and only this file did not. So the
loop body absorbs s_relax_cell_pressure, which existed only to hold three per-cell calls, and the
mixture goes back through the shared operator. The file is 38 lines shorter for it.

nvfortran 23.11 and 24.1 segfaulted in fort2 compiling the HLL solver, chemistry config only.
s_compute_chemistry_average_state took molecular_weights - a parameter array from
m_thermochem - as an actual argument into a declare-target routine. Every other use in MFC reads
that array directly in the kernel. It now takes the per-species gas constants, formed at the call
site. Confirmed as an A/B/A on one machine: reverting only this change brings the ICE back, so it is
this and not the accompanying dead-local removal.

CCE OpenMP and OpenACC both faulted the GPU at runtime on the two reacting-Roe chemistry cases -
Memory access fault by GPU node-N ... Reason: Unknown on the first step, with a clean build.
s_compute_chemistry_average_state called m_thermochem one routine deeper than master did, holding
seven num_species automatic arrays in the intermediate frame. This is the same rule as the
num_fluids arrays above, in a second guise: such calls belong in the parallel-loop body. The
enthalpies and heat capacities are now evaluated at the call sites and passed in, leaving the routine
pure arithmetic over three arrays - and both call sites already computed the heat capacities, so it
had been recomputing them. Chem-config GPU resources move by two registers on one kernel of 473.

Note that cray_inline / cray_noinline / function_name on GPU_ROUTINE are no-ops on the
OpenACC path - all three branches emit the same !$acc routine, and the !DIR$ hints appear only on
Cray CPU builds. Neither of these is fixable by an inlining hint.

Fix: the Roe average dropped the transverse kinetic energy

Found while centralizing the sound speed. roe_avg accumulated vel_avg_rms over all velocity
components and then overwrote it with the first component alone:

do i = 1, num_vels
    vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp
end do
...
vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp

vel_avg_rms has exactly one consumer, s_compute_speed_of_sound_avg, which needs the full squared
magnitude of the averaged velocity. Four things agree the loop is right and the overwrite is wrong:

1. c = (H - 0.5*vel_sum - qv/rho)/gamma is the stiffened-gas Roe sound speed, in which vel_sum is
   |u|^2 summed over all components.
2. It is not even the normal component. dir_idx is (1,2,3), (2,1,3), (3,1,2) per sweep and
   momentum stays in physical order, so vel_L(dir_idx(1)) is the normal velocity and vel_L(1) is
   the x-component in every sweep - a transverse velocity in the y and z sweeps.
3. arithmetic_avg keeps the full sum, and so does the hand-written bubbles path at
   m_riemann_solver_hllc.fpp:625. Only roe_avg overwrote.
4. In the chemistry branch vel_avg_rms appears twice with opposite sign - once inside Phi_avg,
   once as vel_sum - and the two cancel to leave c^2 = gamma*R*T exactly, but only for the full
   magnitude. With the overwrite the residue is + (gamma-1)*|u_transverse|^2/2.

Deleting the overwrite line is the fix. The error is zero in 1D (num_vels == 1), which is why the
two flamelet examples that use avg_state = 1 with wave_speeds = 2 never showed it.

Why it survived: roe_avg runs only when avg_state == 1 and wave_speeds == 2, and
cases.py pushed those as sibling cases off the same stack level, never combined - so the
avg_state=1 cases computed the average and discarded it, and the wave_speeds=2 cases took the
arithmetic branch. The path had no effective coverage. This PR adds the missing combination for 1D,
2D and 3D against HLL and HLLC.

Generating those six goldens on the unfixed code and then regenerating after the fix gives exactly
the predicted split: the two 1D goldens come back byte-identical, all four 2D/3D goldens move.
No pre-existing golden changes - roe_avg is reachable only through compute_average_state, which
is called only under wave_speeds == wave_speeds_pressure, so {avg_state == 1 and wave_speeds == 2} is the complete affected set and nothing else in the tree is in it.

Verification: 97 targeted cases covering every riemann_solver, avg_state, wave_speeds,
low_Mach and mixture_err combination across 1D/2D/3D and both fluid counts - 97 passed, 0
failed. Static GPU resources 0 regression(s) across 471 kernels for the Roe fix and again for the
low-Mach conversion, each measured as a before/after pair on one node.

Fix: the hypoelastic strain energy summed one stress component for all of them (#1773)

s_compute_pressure declared stress as a scalar and then added it once per stress index, so the
elastic energy was 4*tau_11^2 instead of tau_11^2 + 2*tau_12^2 + tau_22^2. Every caller passes
eqn_idx%stress%beg. Correct in 1D, where the one component is not a shear one; wrong in 2D and 3D.

The caller now supplies the summed energy, which removes the loop, shear_indices and a loop
iterator from a device routine in src/common. f_elastic_energy already computed the per-component
term, so it moves from m_riemann_state to m_variables_conversion beside the other operators and
gains callers rather than copies. The conservative-to-primitive conversion carried a fourth copy in
the per-cell GPU path; it now calls the same function, which also absorbs its near-zero shear-modulus
guard and its shear doubling - ten lines become three. Four copies become one, gated on that hot path
by unchanged goldens and 0 regression(s) across 471 kernels.

Only probe output reaches this branch, and no test enabled probe_wrt on a multi-dimensional
hypoelastic case, so the fix brings one: 2D, one fluid, tau_e = (-1e4, 3e3, -2e3). Equal components
would hide the bug - both forms give 4*tau^2 - so the stresses have to differ.

The corrected value is the pressure the case specifies, a 2.4% error removed, and it matches a
hand-evaluation of (E_buggy - E_correct)/gamma to every printed digit. No other field in that
golden moves, and the three 1D hypoelastic cases are unchanged, as the 1D-is-correct argument
predicts.

Fix: HLLC never used the reacting Roe sound speed (#1774)

Two halves. HLL read c_sum_Yi_Phi uninitialized on the arithmetic path - .and. is not required to
short-circuit, so the read happens even though the value is unused. It is now zeroed before the
branch.

The other half is a wrong answer rather than a latent one. The HLLC base path formed the chemistry
average state and then passed a literal 0._wp; since abs(c_c) > verysmall gates the branch, it
always fell back to the frozen sound speed, while HLL and the hypoelastic variant used the Roe one.
The same case got a different sound speed depending on which solver ran it and whether
hypoelasticity was on.

The two are not interchangeable. c_c - (gamma - 1)*(vel_sum - H) reduces exactly to the Roe
sound speed (gamma - 1)*(H_avg - |u_avg|^2/2) - the (gamma - 1)*h_avg and gamma*R*T_avg terms
cancel - whereas the fallback is the frozen mixture speed evaluated at pres_R with rho_avg (which
the original comment called placeholders) and drops the velocity-variance term the linearization
needs. On a representative shocktube state they differ by ~2%. avg_state = 1 asks for Roe
averaging, so the Roe sound speed is what it should get.

The bubbles_euler loop keeps its literal zero and now says why: it never forms the chemistry
average, and chemistry with bubbles_euler/qbmm is prohibited, so that branch is unreachable.
Replacing that zero too would have reintroduced the uninitialized read.

Of the fourteen chemistry cases, only the HLL case added here reached the Roe branch at all, so the
HLLC combination had no coverage and no existing golden moves. A riemann_solver = 2 case is added
beside it.